﻿Imports System
Imports System.IO
Imports System.Runtime.Serialization
Imports System.Runtime.Serialization.Formatters.Binary
''' <summary>
''' dem文件管理类，用于管理、合并、修改、生成DEM文件
''' </summary>
''' <remarks></remarks>
<Serializable()> Public Class DemManage
    ''' <summary>
    ''' DEM文件名
    ''' </summary>
    ''' <remarks></remarks>
    Public FileName As String = ""
    ''' <summary>
    ''' 标准的初始化头文件，可直接在这个基础上进行修改
    ''' </summary>
    ''' <remarks></remarks>
    Public IntialHead As String
    Private m_PlaniMetricreference As Integer  '坐标系统。0表示GEO坐标，1表示UTM坐标
    Private m_Zone As Integer 'UTM区。当坐标系为GEO时，为0。
    Private m_Unit As Integer '单位。坐标系的单位0=radians、1=feet、2=meters、3=arc-seconds
    Private m_ElevationUnit As Integer '高程的单位。1=feet、2=meters
    Private m_CoordinatesSW As New Point3D 'DEM数据西南角的坐标
    Private m_CoordinatesNW As New Point3D 'DEM数据西北角的坐标
    Private m_CoordinatesNE As New Point3D 'DEM数据东北角的坐标
    Private m_CoordinatesSE As New Point3D 'DEM数据东南角的坐标
    Private m_xWidth As Double '数据在X轴方向上的宽度
    Private m_yWidth As Double '数据在Y轴上的宽度
    Private m_ZWidth As Double '数据在Z轴上的系数。即程的系数。比如0.1，即表示用高程数据*0.1后得到高度。
    Private m_Minimum As Double 'DEM数据中的最小值
    Private m_Maximum As Double 'DEM数据中的最大值
    Private m_Rows As Integer 'DEM数据的行。指的是DEM图上的行
    Private m_Cols As Integer 'DEM数据的列。指的是DEM图上的列
    Public m_A(-1, -1) As Integer 'DEM数据区域。按图上的行列，与文本文件是反的
    Public m_ARow(-1) As Integer '用于储存图上每一列上有几行
    ''' <summary>
    ''' 坐标系统。0表示GEO坐标，1表示UTM坐标
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property PlaniMetricreference() As Integer
        Get
            Return m_PlaniMetricreference
        End Get
        Set(ByVal value As Integer)
            Me.m_PlaniMetricreference = value
        End Set
    End Property
    ''' <summary>
    ''' UTM区。当坐标系为GEO时，为0。
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property Zone() As Integer
        Get
            Return m_Zone
        End Get
        Set(ByVal value As Integer)
            m_Zone = value
        End Set
    End Property
    ''' <summary>
    ''' 单位。坐标系的单位0=radians、1=feet、2=meters、3=arc-seconds
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property Unit() As Integer
        Get
            Return m_Unit
        End Get
        Set(ByVal value As Integer)
            m_Unit = value
        End Set
    End Property
    ''' <summary>
    ''' 高程的单位。1=feet、2=meters
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property ElevationUnit() As Integer
        Get
            Return Me.m_ElevationUnit
        End Get
        Set(ByVal value As Integer)
            Me.m_ElevationUnit = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据西南角的坐标
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property CoordinatesSW() As Point3D
        Get
            Return m_CoordinatesSW
        End Get
        Set(ByVal value As Point3D)
            m_CoordinatesSW = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据西北角的坐标
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property CoordinatesNW() As Point3D
        Get
            Return m_CoordinatesNW
        End Get
        Set(ByVal value As Point3D)
            m_CoordinatesNW = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据东北角的坐标
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property CoordinatesNE() As Point3D
        Get
            Return m_CoordinatesNE
        End Get
        Set(ByVal value As Point3D)
            m_CoordinatesNE = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据东南角的坐标
    ''' </summary>
    ''' <remarks></remarks>
    Property CoordinatesSE() As Point3D
        Get
            Return m_CoordinatesSE
        End Get
        Set(ByVal value As Point3D)
            m_CoordinatesSE = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据中的最小值
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property Minimum() As Double
        Get
            Return m_Minimum
        End Get
        Set(ByVal value As Double)
            m_Minimum = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据中的最大值
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property Maximum() As Double
        Get
            Return m_Maximum
        End Get
        Set(ByVal value As Double)
            m_Maximum = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据的行
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property Rows() As Integer
        Get
            Return m_Rows
        End Get
        Set(ByVal value As Integer)
            m_Rows = value
        End Set
    End Property
    ''' <summary>
    ''' DEM数据的列
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property Cols() As Integer
        Get
            Return m_Cols
        End Get
        Set(ByVal value As Integer)
            m_Cols = value
        End Set
    End Property
    ''' <summary>
    ''' 数据在X轴方向上的宽度
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property xWidth() As Double
        Get
            Return m_xWidth
        End Get
        Set(ByVal value As Double)
            m_xWidth = value
        End Set
    End Property
    ''' <summary>
    ''' 数据在Y轴上的宽度
    ''' </summary>
    ''' <value></value>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Property yWidth() As Double
        Get
            Return m_yWidth
        End Get
        Set(ByVal value As Double)
            m_yWidth = value
        End Set
    End Property
    ''' <summary>
    ''' 数据在Z轴上的系数。即程的系数。比如0.1，即表示用高程数据*0.1后得到高度。
    ''' </summary>
    ''' <remarks></remarks>
    Property ZWidth() As Double
        Get
            Return m_ZWidth
        End Get
        Set(ByVal value As Double)
            m_ZWidth = value
        End Set
    End Property
    ''' <summary>
    ''' 获取DEM数据头文件的信息
    ''' </summary>
    ''' <param name="strall"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Private Function GetDemHeadFormString(ByVal strall As String) As Boolean
        Try
            IntialHead = strall.Substring(0, 1024) '取得头文件
            m_PlaniMetricreference = IntialHead.Substring(156, 6)   '坐标系统。0表示GEO坐标，1表示UTM坐标
            m_Zone = IntialHead.Substring(162, 6) 'UTM区。当坐标系为GEO时，为0。
            m_Unit = IntialHead.Substring(528, 6) '单位。坐标系的单位0=radians、1=feet、2=meters、3=arc-seconds
            m_ElevationUnit = IntialHead.Substring(534, 6) '高程的单位。1=feet、2=meters
            '读取4个角的坐标
            Dim strLongSWx As String = IntialHead.Substring(546, 24) '西南角的经度
            Dim strLongSWy As String = IntialHead.Substring(570, 24) '西南角的经度
            Dim strLongNWx As String = IntialHead.Substring(594, 24) '西北角的经度
            Dim strLongNWy As String = IntialHead.Substring(618, 24) '西北角的经度
            Dim strLongNEx As String = IntialHead.Substring(642, 24) '东北角的经度
            Dim strLongNEy As String = IntialHead.Substring(666, 24) '东北角的经度
            Dim strLongSEx As String = IntialHead.Substring(690, 24) '东北角的经度
            Dim strLongSEy As String = IntialHead.Substring(714, 24) '东北角的经度
            strLongSWx = strLongSWx.Replace("D", "E")
            strLongNWx = strLongNWx.Replace("D", "E")
            strLongNEx = strLongNEx.Replace("D", "E")
            strLongSEx = strLongSEx.Replace("D", "E")
            strLongSWy = strLongSWy.Replace("D", "E")
            strLongNWy = strLongNWy.Replace("D", "E")
            strLongNEy = strLongNEy.Replace("D", "E")
            strLongSEy = strLongSEy.Replace("D", "E")
            m_CoordinatesSW.x = CDbl(strLongSWx)
            m_CoordinatesSW.y = CDbl(strLongSWy)
            m_CoordinatesNW.x = CDbl(strLongNWx)
            m_CoordinatesNW.y = CDbl(strLongNWy)
            m_CoordinatesNE.x = CDbl(strLongNEx)
            m_CoordinatesNE.y = CDbl(strLongNEy)
            m_CoordinatesSE.x = CDbl(strLongSEx)
            m_CoordinatesSE.y = CDbl(strLongSEy)
            '读取宽度和高度尺度
            Dim strxWidth As String = IntialHead.Substring(816, 12)
            Dim stryWidth As String = IntialHead.Substring(828, 12)
            Dim strzWidth As String = IntialHead.Substring(840, 12)
            strxWidth = strxWidth.Replace("D", "E")
            stryWidth = stryWidth.Replace("D", "E")
            strzWidth = strzWidth.Replace("D", "E")
            m_xWidth = CDbl(strxWidth) '数据在X轴方向上的宽度
            m_yWidth = CDbl(stryWidth) '数据在Y轴上的宽度
            m_ZWidth = CDbl(strzWidth) '数据在Z轴上的系数。即程的系数。比如0.1，即表示用高程数据*0.1后得到高度。
            Dim strMinimum As String = IntialHead.Substring(738, 24)
            Dim strMaximum As String = IntialHead.Substring(762, 24)
            strMinimum = strMinimum.Replace("D", "E")
            strMaximum = strMaximum.Replace("D", "E")
            m_Minimum = CDbl(strMinimum) 'DEM数据中的最小值
            m_Maximum = CDbl(strMaximum) 'DEM数据中的最大值
            m_Rows = IntialHead.Substring(852, 6) 'DEM数据的行
            m_Cols = IntialHead.Substring(858, 6) 'DEM数据的列
            Return True
        Catch ex As Exception
            Return False
        End Try
    End Function
    Public Function GetDemFormString(ByVal strall As String) As Boolean
        Try
            IntialHead = strall.Substring(0, 1024) '取得头文件
            m_PlaniMetricreference = IntialHead.Substring(156, 6)   '坐标系统。0表示GEO坐标，1表示UTM坐标
            m_Zone = IntialHead.Substring(162, 6) 'UTM区。当坐标系为GEO时，为0。
            m_Unit = IntialHead.Substring(528, 6) '单位。坐标系的单位0=radians、1=feet、2=meters、3=arc-seconds
            m_ElevationUnit = IntialHead.Substring(534, 6) '高程的单位。1=feet、2=meters
            '读取4个角的坐标
            Dim strLongSWx As String = IntialHead.Substring(546, 24) '西南角的经度
            Dim strLongSWy As String = IntialHead.Substring(570, 24) '西南角的经度
            Dim strLongNWx As String = IntialHead.Substring(594, 24) '西北角的经度
            Dim strLongNWy As String = IntialHead.Substring(618, 24) '西北角的经度
            Dim strLongNEx As String = IntialHead.Substring(642, 24) '东北角的经度
            Dim strLongNEy As String = IntialHead.Substring(666, 24) '东北角的经度
            Dim strLongSEx As String = IntialHead.Substring(690, 24) '东北角的经度
            Dim strLongSEy As String = IntialHead.Substring(714, 24) '东北角的经度
            strLongSWx = strLongSWx.Replace("D", "E")
            strLongNWx = strLongNWx.Replace("D", "E")
            strLongNEx = strLongNEx.Replace("D", "E")
            strLongSEx = strLongSEx.Replace("D", "E")
            strLongSWy = strLongSWy.Replace("D", "E")
            strLongNWy = strLongNWy.Replace("D", "E")
            strLongNEy = strLongNEy.Replace("D", "E")
            strLongSEy = strLongSEy.Replace("D", "E")
            m_CoordinatesSW.x = CDbl(strLongSWx)
            m_CoordinatesSW.y = CDbl(strLongSWy)
            m_CoordinatesNW.x = CDbl(strLongNWx)
            m_CoordinatesNW.y = CDbl(strLongNWy)
            m_CoordinatesNE.x = CDbl(strLongNEx)
            m_CoordinatesNE.y = CDbl(strLongNEy)
            m_CoordinatesSE.x = CDbl(strLongSEx)
            m_CoordinatesSE.y = CDbl(strLongSEy)
            '读取宽度和高度尺度
            Dim strxWidth As String = IntialHead.Substring(816, 12)
            Dim stryWidth As String = IntialHead.Substring(828, 12)
            Dim strzWidth As String = IntialHead.Substring(840, 12)
            strxWidth = strxWidth.Replace("D", "E")
            stryWidth = stryWidth.Replace("D", "E")
            strzWidth = strzWidth.Replace("D", "E")
            m_xWidth = CDbl(strxWidth) '数据在X轴方向上的宽度
            m_yWidth = CDbl(stryWidth) '数据在Y轴上的宽度
            m_ZWidth = CDbl(strzWidth) '数据在Z轴上的系数。即程的系数。比如0.1，即表示用高程数据*0.1后得到高度。
            Dim strMinimum As String = IntialHead.Substring(738, 24)
            Dim strMaximum As String = IntialHead.Substring(762, 24)
            strMinimum = strMinimum.Replace("D", "E")
            strMaximum = strMaximum.Replace("D", "E")
            m_Minimum = CDbl(strMinimum) 'DEM数据中的最小值
            m_Maximum = CDbl(strMaximum) 'DEM数据中的最大值
            m_Rows = IntialHead.Substring(852, 6) 'DEM数据的行
            m_Cols = IntialHead.Substring(858, 6) 'DEM数据的列
            '如果Rows为1,cols表示DEM数据有几列
            '可以通过计算求出最多会有几行。
            Dim n As Integer = (m_CoordinatesNE.y - m_CoordinatesSW.y) / yWidth + 1
            Me.Rows = n
            ReDim m_A(n - 1, m_Cols - 1) 'DEM数据区域。也按DEM图上的行、列对应。行是从下向上，列是从左向右
            For i As Integer = 0 To m_Rows - 1
                For j As Integer = 0 To m_Cols - 1
                    m_A(i, j) = -32767
                Next
            Next
            ReDim m_ARow(m_Cols - 1) '用于储存每一列上有几行

            '由于每一个DEM数据的记录由多行结成，在DEM文件中行中储存的是图中列的数据。先判断图中的一列数据储存在几行文本中
            Dim nown As Integer = 1024 '用于计算读到第几行了。每一行1024个字符。第0个已经被作为head行读取。所以从第1025个开始。
            For j As Integer = 0 To m_Cols - 1
                Dim i As Integer = 0
                '先读取第一列有几行
                Dim strLine As String = strall.Substring(nown, 1024)
                m_ARow(j) = strLine.Substring(12, 6)
                '先读取每个记录的第一行
                Dim fisrtrow As Integer
                If m_ARow(j) > 146 Then
                    fisrtrow = 146
                Else
                    fisrtrow = m_ARow(j)
                End If
                For k As Integer = 0 To fisrtrow - 1
                    m_A(i, j) = strLine.Substring(144 + 6 * k, 6)
                    i += 1
                Next
                nown += 1024
                If m_ARow(j) > 146 Then
                    '每6格一个数
                    Dim r As Integer = Math.Truncate((m_ARow(j) - 146 - 1) / 170) + 1
                    For L As Integer = 0 To r - 1
                        strLine = strall.Substring(nown, 1024)
                        For k As Integer = 0 To 170 - 1
                            If strLine.Substring(k * 6, 6) <> "      " Then
                                m_A(i, j) = strLine.Substring(k * 6, 6)
                            End If
                            i += 1
                            If i > m_ARow(j) - 1 Then
                                Exit For
                            End If
                        Next
                        nown += 1024
                    Next
                End If
            Next
            Return True
        Catch ex As Exception
            Return False
        End Try
    End Function


    ''' <summary>
    ''' 读取DEM文件。如果成功返回True，否则返回false
    ''' </summary>
    ''' <param name="sfilename">文件名</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Function ReadDEM(ByVal sFileName As String) As Boolean

        Dim strAll As String = My.Computer.FileSystem.ReadAllText(sFileName)
        Return GetDemFormString(strAll)
    End Function


    ''' <summary>
    ''' 读取DEM文件。如果成功返回True，否则返回false
    ''' </summary>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Function ReadDEMHead() As Boolean
        'Try
        '    If My.Computer.FileSystem.FileExists(Me.ProjPath & "A0.dem") Then
        '        My.Computer.FileSystem.DeleteFile(Me.ProjPath & "A0.dem")
        '    End If
        '    My.Computer.FileSystem.CopyFile(dlgOpen.FileName, Me.ProjPath & "A0.dem")
        'Catch ex As Exception
        '    MsgBox("文件正在被使用")
        '    Exit Function
        'End Try

        Dim strAll As String = My.Computer.FileSystem.ReadAllText(Me.FileName)
        Return GetDemHeadFormString(strAll)
    End Function
    ''' <summary>
    ''' 写DEM文件。如果成功返回True，否则返回false
    ''' </summary>
    ''' <param name="sfilename"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Function WriteDEM(ByVal sFileName As String) As Boolean
        '根据用户的设置修改头文件

        IntialHead = IntialHead.Remove(162, 6)
        IntialHead = IntialHead.Insert(162, m_Zone.ToString.PadLeft(6))

        IntialHead = IntialHead.Remove(546, 24) '西南角的经度
        IntialHead = IntialHead.Insert(546, m_CoordinatesSW.x.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(570, 24) '西南角的经度
        IntialHead = IntialHead.Insert(570, m_CoordinatesSW.y.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(594, 24) '西北角的经度
        IntialHead = IntialHead.Insert(594, m_CoordinatesNW.x.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(618, 24) '西北角的经度
        IntialHead = IntialHead.Insert(618, m_CoordinatesNW.y.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(642, 24) '东北角的经度
        IntialHead = IntialHead.Insert(642, m_CoordinatesNE.x.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(666, 24) '东北角的经度
        IntialHead = IntialHead.Insert(666, m_CoordinatesNE.y.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(690, 24) '东北角的经度
        IntialHead = IntialHead.Insert(690, m_CoordinatesSE.x.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(714, 24) '东北角的经度
        IntialHead = IntialHead.Insert(714, m_CoordinatesSE.y.ToString.PadLeft(24))


        IntialHead = IntialHead.Remove(738, 24)
        IntialHead = IntialHead.Insert(738, m_Minimum.ToString.PadLeft(24))
        IntialHead = IntialHead.Remove(762, 24)
        IntialHead = IntialHead.Insert(762, m_Maximum.ToString.PadLeft(24))

        IntialHead = IntialHead.Remove(852, 6)
        IntialHead = IntialHead.Insert(852, m_Rows.ToString.PadLeft(6))
        IntialHead = IntialHead.Remove(858, 6)
        IntialHead = IntialHead.Insert(858, m_Cols.ToString.PadLeft(6))

        Try
            Dim strAll As String = IntialHead
            For j As Integer = 0 To m_Cols - 1
                Dim i As Integer = 0
                '先找出每一列的最大值和最小值
                Dim min As Double
                Dim max As Double
                Dim getmaxandmin As Boolean = False
                For l As Integer = 0 To m_ARow(j) - 1
                    If m_A(l, j) < -30000 Then
                        min = m_A(l, j)
                        max = m_A(l, j)
                        getmaxandmin = True
                        Exit For
                    End If
                Next
                If getmaxandmin = True Then
                    For l As Integer = 0 To m_ARow(j) - 1
                        If m_A(l, j) > max Then
                            max = m_A(l, j)
                        End If
                        If m_A(l, j) < min Then
                            min = m_A(l, j)
                        End If
                    Next
                End If
                '先写第一列
                Dim strLine As String = "     1" & (j + 1).ToString.PadLeft(6) & m_ARow(j).ToString.PadLeft(6) & "     1" _
                                        & (m_CoordinatesSW.x + m_xWidth * j).ToString.PadLeft(24) & m_CoordinatesSW.y.ToString.PadLeft(24) _
                                        & "0.000000000000000".PadLeft(24) & min.ToString.PadLeft(24) & max.ToString.PadLeft(24)
                Dim fisrtrow As Integer
                If m_ARow(j) > 146 Then
                    fisrtrow = 146
                Else
                    fisrtrow = m_ARow(j)
                End If
                For k As Integer = 0 To fisrtrow - 1
                    strLine = strLine & m_A(i, j).ToString.PadLeft(6)
                    i += 1
                Next
                strLine = strLine.PadRight(1024)
                If m_ARow(j) > 146 Then
                    '每6格一个数
                    Dim r As Integer = Math.Truncate((m_ARow(j) - 146 - 1) / 170) + 1
                    For L As Integer = 0 To r - 1
                        For k As Integer = 0 To 170 - 1
                            strLine = strLine & m_A(i, j).ToString.PadLeft(6)
                            i += 1
                            If i > m_ARow(j) - 1 Then
                                Exit For
                            End If
                        Next
                        strLine = strLine.PadRight(1024 * (L + 2))
                    Next
                End If
                strAll = strAll & strLine
            Next
            'strAll = strAll.Replace(" ", Chr(32))
            'My.Computer.FileSystem.WriteAllText(sFileName, strAll, False)
            Dim abyte(strAll.Length - 1) As Byte
            For j As Integer = 0 To strAll.Length - 1
                abyte(j) = AscW(strAll(j))
            Next
            My.Computer.FileSystem.WriteAllBytes(sFileName, abyte, False)
        Catch ex As Exception
            Return False
        End Try
        Return True
    End Function

    Public Shared Function UniteDem(ByVal dem1 As DemManage, ByVal dem2 As DemManage) As DemManage
        '先判断两个文件的位置
        '1    2
        '4    3
        '先找出两个文件中的最大和最小坐标
        Dim minx, miny, maxx, maxy As Double
        minx = dem1.CoordinatesSW.x
        miny = dem1.CoordinatesSW.y
        maxx = dem1.CoordinatesNE.x
        maxy = dem1.m_CoordinatesNE.y
        If dem2.CoordinatesSW.x < minx Then
            minx = dem2.CoordinatesSW.x
        End If
        If dem2.CoordinatesSW.y < miny Then
            miny = dem2.CoordinatesSW.y
        End If
        If dem2.CoordinatesNE.x > maxx Then
            maxx = dem2.CoordinatesNE.x
        End If
        If dem2.CoordinatesNE.y > maxy Then
            maxy = dem2.CoordinatesNE.y
        End If
        Dim NewDem As New DemManage



        NewDem.IntialHead = dem1.IntialHead '取得头文件
        '读取4个角的坐标
        'NewDem.CoordinatesSW.x = minx
        'NewDem.CoordinatesSW.y = miny
        'NewDem.CoordinatesNW.x = minx
        'NewDem.CoordinatesNW.y = maxy
        'NewDem.CoordinatesNE.x = maxx
        'NewDem.CoordinatesNE.y = maxy
        'NewDem.CoordinatesSE.x = maxx
        'NewDem.CoordinatesSE.y = miny
        '最大值和最小值
        If dem1.Minimum < dem2.Minimum Then
            NewDem.Minimum = dem1.Minimum
        Else
            NewDem.Minimum = dem2.Minimum
        End If
        If dem1.Maximum > dem2.Maximum Then
            NewDem.Maximum = dem1.Maximum
        Else
            NewDem.Maximum = dem2.Maximum
        End If
        NewDem.yWidth = dem1.yWidth
        NewDem.xWidth = dem1.xWidth
        NewDem.Rows = (maxy - miny) / NewDem.yWidth + 1 'DEM数据的行
        NewDem.Cols = (maxx - minx) / NewDem.xWidth + 1 'DEM数据的列

        '如果Rows为1,cols表示DEM数据有几列
        '可以通过计算求出最多会有几行。
        Dim n As Integer = (NewDem.CoordinatesNE.x - NewDem.CoordinatesSW.x) / dem1.xWidth + 1
        Dim m As Integer = (NewDem.CoordinatesNE.y - NewDem.CoordinatesSW.y) / dem1.yWidth + 1
        ReDim NewDem.m_A(m - 1, n - 1) 'DEM数据区域。也按DEM图上的行、列对应。行是从下向上，列是从左向右
        For i As Integer = 0 To m - 1
            For j As Integer = 0 To n - 1
                NewDem.m_A(i, j) = -32767
            Next
        Next
        ReDim NewDem.m_ARow(n - 1) '用于储存每一列上有几行
        For j As Integer = 0 To n - 1
            NewDem.m_ARow(j) = m
        Next
        '由于每一个DEM数据的记录由多行结成，在DEM文件中行中储存的是图中列的数据。先判断图中的一列数据储存在几行文本中
        Dim mstart1 As Integer = (dem1.CoordinatesSW.y - NewDem.CoordinatesSW.y) / dem1.yWidth
        Dim nstart1 As Integer = (dem1.CoordinatesSW.x - NewDem.CoordinatesSW.x) / dem1.xWidth

        For j As Integer = 0 To dem1.Cols - 1
            For i As Integer = 0 To dem1.m_ARow(j) - 1
                NewDem.m_A(mstart1 + i, nstart1 + j) = dem1.m_A(i, j)
            Next
        Next
        Dim mstart2 As Integer = (dem2.CoordinatesSW.y - NewDem.CoordinatesSW.y) / dem2.yWidth
        Dim nstart2 As Integer = (dem2.CoordinatesSW.x - NewDem.CoordinatesSW.x) / dem2.xWidth
        For j As Integer = 0 To dem2.Cols - 1
            For i As Integer = 0 To dem2.m_ARow(j) - 1
                NewDem.m_A(mstart2 + i, nstart2 + j) = dem2.m_A(i, j)
            Next
        Next
        '判断中间是否有间隙
        If dem1.Cols + dem2.Cols < NewDem.Cols Then '列有间隙
            '间隙的位置在nstart2-1处
            Dim jjj As Integer
            If nstart2 > nstart1 Then
                jjj = nstart2 - 1
            Else
                jjj = nstart1 - 1
            End If
            For i As Integer = 0 To NewDem.Rows - 1
                If NewDem.m_A(i, jjj - 1) > -30000 And NewDem.m_A(i, jjj + 1) > -30000 Then
                    NewDem.m_A(i, jjj) = (NewDem.m_A(i, jjj - 1) + NewDem.m_A(i, jjj + 1)) / 2
                End If
            Next

        End If

        If dem1.Rows + dem2.Rows < NewDem.Rows Then '行有间隙
            '间隙的位置在nstart2-1处
            Dim iii As Integer
            If mstart2 > mstart1 Then
                iii = mstart2 - 1
            Else
                iii = mstart1 - 1
            End If
            For j As Integer = 0 To NewDem.Cols - 1
                If NewDem.m_A(iii - 1, j) > -30000 And NewDem.m_A(iii + 1, j) > -30000 Then
                    NewDem.m_A(iii, j) = (NewDem.m_A(iii - 1, j) + NewDem.m_A(iii + 1, j)) / 2
                End If
            Next
        End If

        Return NewDem
    End Function




    ''' <summary>
    ''' 获得任意一点的Z值
    ''' </summary>
    ''' <remarks></remarks>
    Public Function GetZ(ByVal X As Double, ByVal Y As Double) As Double
        Dim DemPt(15) As Point3D  'pt的16个邻点
        For k As Integer = 0 To 15
            DemPt(k) = New Point3D
        Next


        Dim StartX, StartY As Double
        StartX = Me.m_CoordinatesSW.x
        StartY = Me.m_CoordinatesSW.y

        Dim i, j As Integer

        Dim M_Z As Double

        If X < Me.m_CoordinatesSW.x Or X > Me.m_CoordinatesNE.x Or Y < StartY Or Y > Me.m_CoordinatesNE.y Then
            M_Z = -99
            'MsgBox("坐标超出地形范围，请重新设置！")
            'Exit Function

        ElseIf (X < Me.m_CoordinatesSW.x + Me.xWidth Or X > Me.m_CoordinatesNE.x - Me.xWidth Or Y < StartY + Me.yWidth Or Y > Me.m_CoordinatesNE.y - Me.yWidth) _
         And (X > Me.m_CoordinatesSW.x Or X < Me.m_CoordinatesNE.x Or Y > StartY Or Y < Me.m_CoordinatesNE.y) Then '在边界区域用4个点插值,线性插值
            For i = 0 To Me.m_Cols - 1
                If StartX > X Then

                    DemPt(0).x = StartX - Me.m_xWidth
                    DemPt(1).x = StartX - Me.m_xWidth
                    DemPt(2).x = StartX
                    DemPt(3).x = StartX

                    Exit For
                End If
                StartX = StartX + Me.m_xWidth
            Next


            For j = 0 To Me.m_Rows - 1
                If StartY > Y Then

                    DemPt(0).y = StartY - Me.m_yWidth
                    DemPt(1).y = StartY
                    DemPt(2).y = StartY
                    DemPt(3).y = StartY - Me.m_yWidth

                    Exit For
                End If
                StartY = StartY + Me.m_yWidth
            Next

            DemPt(0).z = Me.m_A(j - 1, i - 1) * Me.m_ZWidth
            DemPt(1).z = Me.m_A(j, i - 1) * Me.m_ZWidth
            DemPt(2).z = Me.m_A(j, i) * Me.m_ZWidth
            DemPt(3).z = Me.m_A(j - 1, i) * Me.m_ZWidth

            M_Z = insert1(X, Y, DemPt)

        Else   '用16个点插值
            For i = 0 To Me.m_Cols - 1
                If StartX > X Then

                    DemPt(0).x = StartX - Me.m_xWidth * 2
                    DemPt(1).x = StartX - Me.m_xWidth * 2
                    DemPt(2).x = StartX - Me.m_xWidth * 2
                    DemPt(3).x = StartX - Me.m_xWidth * 2

                    DemPt(4).x = StartX - Me.m_xWidth
                    DemPt(5).x = StartX - Me.m_xWidth
                    DemPt(6).x = StartX - Me.m_xWidth
                    DemPt(7).x = StartX - Me.m_xWidth

                    DemPt(8).x = StartX
                    DemPt(9).x = StartX
                    DemPt(10).x = StartX
                    DemPt(11).x = StartX

                    DemPt(12).x = StartX + Me.m_xWidth
                    DemPt(13).x = StartX + Me.m_xWidth
                    DemPt(14).x = StartX + Me.m_xWidth
                    DemPt(15).x = StartX + Me.m_xWidth

                    Exit For
                End If
                StartX = StartX + Me.m_xWidth
            Next

            For j = 0 To Me.m_Rows - 1
                If StartY > Y Then

                    'DemPt(0).Y = StartY - Me.m_yWidth
                    'DemPt(1).Y = StartY
                    'DemPt(2).Y = StartY
                    'DemPt(3).Y = StartY - Me.m_yWidth
                    DemPt(0).y = StartY - Me.m_yWidth * 2
                    DemPt(1).y = StartY - Me.m_yWidth
                    DemPt(2).y = StartY
                    DemPt(3).y = StartY + Me.m_yWidth

                    DemPt(4).y = StartY + Me.m_yWidth
                    DemPt(5).y = StartY
                    DemPt(6).y = StartY - Me.m_yWidth
                    DemPt(7).y = StartY - Me.m_yWidth * 2

                    DemPt(8).y = StartY - Me.m_yWidth * 2
                    DemPt(9).y = StartY - Me.m_yWidth
                    DemPt(10).y = StartY
                    DemPt(11).y = StartY + Me.m_yWidth

                    DemPt(12).y = StartY + Me.m_yWidth
                    DemPt(13).y = StartY
                    DemPt(14).y = StartY - Me.m_yWidth
                    DemPt(15).y = StartY - Me.m_yWidth * 2


                    Exit For
                End If
                StartY = StartY + Me.m_yWidth
            Next

            '16个点Z值
            ' Z值的二维组好像是列在前, 行在后


            DemPt(0).z = Me.m_A(j - 2, i - 2) * Me.m_ZWidth
            DemPt(1).z = Me.m_A(j - 1, i - 2) * Me.m_ZWidth
            DemPt(2).z = Me.m_A(j, i - 2) * Me.m_ZWidth
            DemPt(3).z = Me.m_A(j + 1, i - 2) * Me.m_ZWidth

            DemPt(7).z = Me.m_A(j - 2, i - 1) * Me.m_ZWidth
            DemPt(6).z = Me.m_A(j - 1, i - 1) * Me.m_ZWidth
            DemPt(5).z = Me.m_A(j, i - 1) * Me.m_ZWidth
            DemPt(4).z = Me.m_A(j + 1, i - 1) * Me.m_ZWidth

            DemPt(8).z = Me.m_A(j - 2, i) * Me.m_ZWidth
            DemPt(9).z = Me.m_A(j - 1, i) * Me.m_ZWidth
            DemPt(10).z = Me.m_A(j, i) * Me.m_ZWidth
            DemPt(11).z = Me.m_A(j + 1, i) * Me.m_ZWidth

            DemPt(15).z = Me.m_A(j - 2, i + 1) * Me.m_ZWidth
            DemPt(14).z = Me.m_A(j - 1, i + 1) * Me.m_ZWidth
            DemPt(13).z = Me.m_A(j, i + 1) * Me.m_ZWidth
            DemPt(12).z = Me.m_A(j + 1, i + 1) * Me.m_ZWidth


            '下面就是根据16个邻居点进行插值求目标点的Z了

            'Dim Z As Double = DemPt(3).Z * (1 - Pt.X) * (1 - Pt.Y) + DemPt(2).Z * (1 - Pt.Y) * Pt.X + DemPt(1).Z * Pt.X * Pt.Y + DemPt(0).Z * (1 - Pt.X) * Pt.Y

            Dim p As Integer = 0
            Dim z_xy(3, 3) As Double
            For k As Integer = 0 To 3
                For l As Integer = 0 To 3
                    z_xy(k, l) = DemPt(p).z
                    p = p + 1
                Next

            Next

            Dim dx As Double = 1          ' Me.m_xWidth - StartX + X
            Dim dy As Double = 1        ' Me.m_yWidth - StartY + Y

            M_Z = insert(dx, dy, z_xy)

        End If

        Return M_Z
    End Function


    ''' <summary>
    ''' 双三次样条曲线插值方法
    ''' </summary>
    ''' <param name="x">数据在X上的宽度</param>
    ''' <param name="y">数据在Y上的宽度</param>
    ''' <param name="z_xy"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Private Function insert(ByVal x As Double, ByVal y As Double, ByVal z_xy(,) As Double) As Double
        Dim a0, a2, a3, b1, b2, b3, c(3) As Double
        For i As Integer = 0 To 3
            a0 = z_xy(0, i) - z_xy(1, i)
            a2 = z_xy(2, i) - z_xy(1, i)
            a3 = z_xy(3, i) - z_xy(1, i)
            b1 = -a0 / 3.0 + a2 - a3 / 6.0           '//求解系数
            b2 = a0 / 2.0 + a2 / 2.0
            b3 = -a0 / 6.0 - a2 / 2.0 + a3 / 6.0
            c(i) = z_xy(1, i) + b1 * x + b2 * x * x + b3 * x * x * x  ' //记录系数
        Next
        a0 = c(0) - c(1)
        a2 = c(2) - c(1)
        a3 = c(3) - c(1)
        b1 = -a0 / 3.0 + a2 - a3 / 6.0  '                //求解系数
        b2 = a0 / 2.0 + a2 / 2.0
        b3 = -a0 / 6.0 - a2 / 2.0 + a3 / 6.0
        Return (c(1) + b1 * y + b2 * y * y + b3 * y * y * y)

    End Function

    ''' <summary>
    ''' 线性插值
    ''' </summary>
    ''' <param name="x"></param>
    ''' <param name="y"></param>
    ''' <param name="pt"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Private Function insert1(ByVal x As Double, ByVal y As Double, ByVal pt() As Point3D) As Double
        Dim z1, z2, z3, z4 As Double
        Dim S As Double = pt(1).Y - pt(0).Y
        Dim SZ As Double = pt(1).Z - pt(0).Z

        z1 = pt(0).Z + (y - pt(0).Y) * SZ / S

        SZ = pt(2).Z - pt(3).Z
        z2 = pt(3).Z + (y - pt(3).Y) * SZ / S

        S = pt(2).X - pt(1).X
        SZ = pt(2).Z - pt(1).Z
        z3 = pt(1).Z + (x - pt(1).X) * SZ / S


        SZ = pt(3).Z - pt(0).Z
        z4 = pt(0).Z + (x - pt(0).X) * SZ / S

        Dim Z As Double = (z1 + z2 + z3 + z4) / 4
        Return Z


    End Function





End Class
